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SECTION  1 
INTRODUCTION 


1.1  BACKGROUND. 

SHIST  (Seismic  Hard  rock  In-situ  Source  Test)  was  a  proposed  prototype  High 
Explosive  (HE)  seismic  source  function  experiment  that  was  designed  to  provide  a 
benchmark  data  set  for  seismic  source  function  models  in  hard  rock  geology.  It  was 
also  intended  to  provide  a  threshold  of  detection  calibration  point  for  regional  seismic 
monitoring.  The  experiment  was  to  consist  of  20  tons  of  a  fully  contained  HE  charge 
(QM-100/R)  detonated  in  a  granite  test  bed.  The  test  cavity  was  to  be  located  at  a 
depth  of  35  meters  and  excavated  into  a  sphere  by  drilling  an  blasting,  with  the  intent  to 
reduce  any  possible  edge  effects  that  could  result  from  the  use  of  a  cylindrical  charge. 
Instrumentation  was  to  be  placed  in  the  hydrodynamic  region  and  in  the  transition  region 
between  the  hydrodynamic  and  seismic  region. 


Concrete  Filled 
Emplacement  Hole 


20  Ton  QM-100/R 
2.9m  Diameter 
Excavated  Sphere 


Over  Drill  Hole 


Ground  Surface 


Weathered  Granite  with 
1.550  m/s  Seismic  Velocity 


< —  1.1  m  Diameter 


-  15  m  Layer 


Jointed  Granite  with 
3,400  m/s  Seismic  Velocity 


Figure  1  -1 .  Proposed  SHIST  test  configuration  and  emplacement  hole  profile. 


1 


The  proposed  SHIST  test  site  is  located  in  Lincoln  County,  New  Mexico,  adjacent 
to  the  Oscura  mountains  in  the  northeast  portion  of  the  White  Sands  Missile  Range 
(WSMR).  The  site  is  in  the  Oscura  Gap  area,  in  the  NW  1/4  of  sec.  16,  T.9  S,  R.6  E, 
and  includes  the  SPERRE  site  where  WES  conducted  HE  cratering  experiments  in 
1979.  The  SHIST  working  point  was  to  be  situated  in  a  region  where  a  combination  of 
age  and  tectogenetic  factors  have  produced  a  granitic  test  bed  that  is  jointed  and 
weathered.  The  site  geology  was  simplified  in  preliminary  scoping  calculations 
presented  here  to  consist  of  a  top  layer  of  highly  weathered  granite  overlaying  a  bed  of 
jointed  and  fractured  granite.  Air-filled  void  (AFV)  porosities  range  from  1-2%  with  In- 
situ  seismic  velocities  ranging  from  1 ,500  to  3,500  m/s  (competent  granite  has  a  seismic 
velocity  on  the  order  of  5000  m/s).  Since  the  proposed  site  is  at  the  foot  of  the  Oscura 
mountains,  the  ground  surface  has  a  slight  slope  to  it  with  a  potential  vertical  variation  of 
25  ft  along  any  one  gage  radial.  The  proposed  site  also  has  known  fault  lines  in  the 
vicinity  of  the  test  bed,  however  the  closest  fault  was  considered  to  be  sufficiently  far 
away  to  be  neglected  in  this  study. 

Currently,  the  SHIST  test  has  been  postponed,  however  Titan  Research  and 
Technology  (TRT)  had  been  tasked  to  develop  preliminary  material  models  for  the  site 
and  conduct  1 D  and  2D  finite-difference  scoping  calculations  to  serve  as  a  starting  point 
should  the  test  (or  other  HE  test  requiring  a  hard  rock)  be  revived  in  the  next  year  or 
two.  The  intent  of  this  study  was  to  investigate  the  details  of  the  shock  wave 
propagation  into  the  granite  site  and  to  determine  the  seismic  transition  point  (the  range 
at  which  the  stress  wave  developed  by  the  HE  source  begins  to  decay  as  an  elastic 
wave)  to  assist  in  future  gage  placement.  Three  dimensional  effects,  such  as  the 
surface  slope  variation  and  ground  fault  lines,  were  neglected  at  this  time.  Should  the 
program  be  reactivated  in  the  future,  these  effects  should  be  examined  in  closer  detail. 

TRT  was  also  tasked  to  install  the  suite  of  CRALE  codes  (ID,  2D)  and  post¬ 
processing  tools  onto  the  DNAFC  Silicon  Graphics  workstation.  The  SHIST  scoping 
study  was  then  to  be  used  to  verify  the  operation  of  these  codes  on  this  platform. 
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1.2  APPROACH. 


The  SHIST  scoping  study  was  completed  as  three  separate  tasks.  Task  1 : 
Evaluate  the  material  properties  of  the  SHIST  site  and  develop  simple  scoping  Equation 
of  State  (EOS)  models  for  jointed  granite.  Due  to  the  uncertainty  in  the  material 
properties  data  for  the  SHIST  site,  the  study  was  extended  to  include  material  model 
excursions  from  the  prescribed  SHIST  lithology.  Task  2:  Take  the  granite  models 
developed  in  Task  1  and  use  them  in  a  series  of  1-D  finite-difference  calculations  (using 
the  CRALE  1 D  code)  to  investigate  the  effect  of  changes  to  material  properties  on  the 
free-field  environment  generated  in  the  testbed.  The  intent  of  these  calculations  was  to 
evaluate  the  potential  effects  material  property  variations  may  have  on  the  seismic 
transition  point.  Task  3:  Take  the  “best-estimate”  material  properties  for  the  site  and 
use  them  in  a  2-D  finite-difference  calculation  (using  the  CRALE  2D  code)  to  investigate 
the  details  of  the  interaction  of  the  shock  with  the  surrounding  stemming  materials  and 
granite.  An  assessment  of  the  degree  of  sphericity  of  the  ground  shock  could  also  be 
made  by  comparing  the  calculated  2-D  results  from  Task  3  with  the  1-D  calculations  of 
Task  2. 

This  report  details  the  modeling  and  calculational  efforts  performed  by  TRT  In 
support  of  the  SHIST  test.  The  material  modeling  is  discussed  Section  2,  the  1-D  free 
field  environment  calculations  in  Section  3,  and  the  2-D  SHIST  scoping  calculations  in 
Section  4.  Appendices  are  included  to  document  the  specifics  of  the  material  models 
used  in  this  study. 

Though  the  proposed  SHIST  experiment  is  fully  coupled,  the  same  analytical 
approach  reported  here  could  be  used  to  characterize  the  wave  in  the  seismic  region  for 
decoupled  sources.  Both  CRALE  1 D  and  2D  codes  could  be  used  to  calculate  the 
shock  wave  resulting  from  an  explosive  source  placed  inside  an  open  cavity.  Surface 
effects  and  reflections  from  layering  could  also  be  examined  In  either  1 D  or  2D 
geometries. 
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SECTION  2 


MATERIAL  MODELS  FOR  GRANITE  AND  CONCRETE 

(TASK1) 

2.1  GENERIC  GRANITE  MODEL  DESCRIPTIONS. 

Prior  to  the  suspension  of  the  SHIST  program,  TRT  developed  a  generic  1% 
(AFV)  weak  granite  equation  of  state  (EOS)  which  was  considered  to  be  the  best  one  to 
use  for  SHIST  p’re-shot  calculations  until  additional  experimental  data  became 
available  [1].  This  model  was  developed  using  an  older  version  of  the  SHEL  EOS  that 
is  no  longer  compatible  with  the  current  version  of  the  CRALE  code  (CRALEPL3). 
Therefore,  under  the  current  effort,  new  granite  material  models  with  characteristics 
similar  to  those  reported  in  Reference  1  were  developed  using  the  current  version  of  the 
SHEL  EOS. 

The  generic  granite  model  discussed  in  Reference  1  (material  model  GRAN3A, 
which  in  this  report  will  be  referred  to  as  GENERIC1)  was  shown  to  produce  a 
reasonable  overall  match  to  data  obtained  from  underground  unclear  tests  at  the 
Nevada  Test  Site  (NTS)  [2].  The  model,  having  a  seismic  velocity  of  4,500  m/s.  Is  one 
representative  of  weak,  jointed  granitic  rock,  and  produces  waveforms  which  agree  well 
with  those  measured  in  the  HARDHAT  test  [2],  peaks  and  scaled  nuclear  cavity  radii 
consistent  with  the  U.S.  Granite  experience  at  NTS  (Rc  =  11 .5  m/kt^'®).  The  generic 
model  incorporates  1%  air  filled  porosity  to  simulate  the  effect  of  natural  joints  and 
fractures  in  otherwise  competent  Granite.  The  failure  surface  parameter  for  GENERIC1 
is  a  simple  linear  Drucker-Prager  /  Von-Mises  yield  surface  of  the  form; 


^  The  current  SHEL  EOS  formulation  is  discussed  in  Appendix  A. 
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(2.1) 


Y  =  max 


'  vm 


where 

=  the  limiting  yield  strength 

Yo  =  the  yield  strength  at  zero  pressure  (cohesion) 
P  =  the  current  pressure,  and 
1/Py  =  the  Drucker-Prager  slope. 


To  maintain  compatibility  with  the  current  CRALE  code,  two  new  generic  granite 
models  (material  models  GRAN1C  and  GRAN1 E,  which  in  this  report  will  be  referred  to 
as  GENERIC2  and  GENERIC3)  were  developed  using  the  current  EOS  model  to  match 
the  characteristics  of  the  earlier  model  (GENERIC1).  Both  of  the  new  models  have  1% 
air  filled  porosity  and  the  identical  maximum  Von-Mises  yield  surface  as  GENERIC1, 
however  two  different  yield  surface  functions  were  used.  The  GENERIC2  uses  the 
same  linear  Drucker-Prager  yield  surface  function  as  GENERIC1  (Equation  2.1). 
GENERIC3  on  the  other  hand  uses  an  exponential  yield  surface  function  of  the  form: 


Y  =  Y,^-(Y 


vm 


■Yo)e‘^^y 


(2.2) 


where 

Y^^  =  the  limiting  yield  strength 

Y^  =  the  yield  strength  at  zero  pressure  (cohesion) 

P  =  the  current  pressure,  and 

Py  =  the  pressure  at  which  the  yield  function  ~(2/3  Y^  +  1/3  Yq)' 

The  exponential  form  of  the  yield  surface  eliminates  the  “non-physical”  sharp 
corner  that  appears  in  Equation  2.1  at  the  intersection  of  the  Drucker-Prager  slope  with 
the  Von-Mises  limit.  Both  formulations  are  available  in  the  current  CRALE  code. 
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The  bulk  properties  and  failure  surface  parameters  for  the  generic  granite  models 
in  Table  2-1  (specific  CRALE  EOS  parameters  for  GENERIC2  and  GENERICS  can  be 
found  in  Appendix  A  and  Appendix  B). 


Table  2-1 .  Properties  of  the  generic  models  for  weak  granite  with  1  %  AFV. 


Model 

Density 

g/cc 

%AFV 

Poisson’s 

ratio 

Cohesion* 

(bar) 

Slope 

1/Py 

Exponential 
Constant  1/Py 

ISiSSSHIi 

2.643 

2.67 

1.0 

0.25 

100 

EH 

N/A 

12.5 

2.624 

2.65 

1.0 

0.25 

100 

0.3 

N/A 

12.5 

Icldjid;ir»kM 

ISIQStnHl 

2.624 

2.65 

1.0 

0.25 

100 

N/A 

30.0 

12.5 

•iPF 


Both  GENERIC2  and  GENERICS  have  a  grain  density  of  2.65  g/cc  (lower  than 
that  of  GENERIC1)  to  be  more  consistent  with  the  SHIST  materiai  properties  data. 

The  uniaxiai  stress-strain  response  of  materials  GENERiCI,  GENERIC2,  and 
GENERICS  are  shown  in  Figures  2-1  a, b  and  2-2a,b,c  below.  The  Hugoniots  for  the 
three  models  are  shown  in  Figure  2-S  along  with  the  currently  available  granite  shock 
data. 
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Figure  2-1 


O 


(a)  stress-strain 


Uniaxial  stress  response  of  GENER1C1  (GRAN3A),  GENERIC2 
(GRAN1C),  and  GENERIC3  (GRAN1E)  to  80  kbar. 
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Figure  2-1 .  Uniaxial  stress  response  from  uniaxial  loading  of  GENERIC1  (GRAN3A), 

GENERIC2  (GRAN1C),  and  GENERIC3  (GRAN1E)  to  80  kbar.  (Continued) 


o 

o 


(a)  stress-strain 


Figure  2-2.  Uniaxial  stress  response  of  GENERIC1  (GRAN3A),  GENERIC2 
(GRAN1C),  and  GENERIC3  (GRAN1E)  to  8  kbar. 
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Figure  2-2.  Uniaxial  stress  response  from  uniaxial  loading  of  GENERIC1  (GRAN3A), 

GENERIC2  (GRAN  1C),  and  GENERIC3  (GRAN IE)  to  80  kbar.  (Continued) 


Figure  2-2  Uniaxial  Stress  response  from  uniaxial  loading  of  GENERIC1  {GRAN3A), 
GENERIC2(GRAN1C),  and  GENERICS  (GRAN1E).  (Continued) 
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Figure  2-3.  Hugoniots  of  GENERIC1  (GRAN3A),  GENERIC2  (GRAN1C), 
and  GENERIC3  (GRAN1 E)  and  representative  experimental 
data  points  to  225  kbar. 
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2.2  GENERIC  GRANITE  MODEL  CHECKOUT  AGAINST  EXISTING 
NUCLEAR  DATA. 


The  peak  velocity  vs.  range  obtained  from  a  fully  contained  1  kt  nuclear 
calculation  employing  the  model  GENERIC1  (GRAN3A)  and  the  two  new  models 
GENERIC2  (GRAN1C)  and  GENERIC3  (GRAN1E)  are  compared  to  available  granite 
test  data  in  Figure  2-4. 


Figure  2-4.  Calculated  peak  velocity  vs.  range  for  GENERIC1  (GRAN3A), 
GENERIC2  (GRAN1C),  and  GENERIC3  (GRAN1E)  models 
compared  to  granite  data  base  (Nuclear  source). 
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As  can  be  seen,  both  new  generic  granite  models  give  virtually  the  same  results 
as  the  older  GENERIC1  model  discussed  in  Reference  1 .  Therefore  these  models,  both 
with  a  seismic  velocity  of  4,500  m/s,  will  be  used  as  material  model  excursions  for  the 
current  SHIST  scoping  study. 

2.3  GRANITE  MATERIAL  MODELS  DEVELOPED  SPECIFICALLY 

FOR  THE  SHIST  SITE. 

The  SHIST  site  is  situated  on  and  adjacent  to  the  pediment  along  the  western 
base  of  the  southern  Oscura  Mountains.  The  Oscura  Mountains  are  an  east-dipping 
north-south  structural  block  with  relatively  gentle  slopes  on  the  east  and  a  steep 
precipitous  scarp  along  the  west.  The  scarp  is  the  result  of  normal  faulting,  an 
expression  of  which  has  been  observed  to  bisect  the  SHIST  site.  For  this  scoping 
study,  the  fault  was  considered  to  be  far  enough  away  from  the  test  bed  so  that  it  could 
be  ignored  at  this  time. 

The  site  profile  in  the  vicinity  of  the  test  bed  is  one  of  highly  weathered  and 
jointed  Precambrian  granite.  Small  sample  laboratory  data  reports  Indicate  that  the  total 
porosities  range  from  1.2  to  1.9%,  with  seismic  velocities  ranging  from  4,300  to  5,300 
m/s  [4].  However,  these  small  scale  properties  are  not  considered  representative  of  the 
in-situ  characteristics  of  the  SHIST  test  bed  since  geophone  data  suggest  that  the 
SHIST  site  has  seismic  velocities  significantly  lower  than  the  lab  data  suggests. 

The  GENERIC  models  discussed  above  have  seismic  velocities  consistent  with 
small  scale  tests  (4,500  m/s)  and  can  be  considered  material  model  excursions  for  the 
SHIST  site  The  actual  SHIST  in-situ  seismic  velocities  range  from  1,200  m/s  in  the 
surface  layers  to  4,000  m/s  in  the  deeper  layers.  Through  previous  discussions  with  Dr. 
Robert  Reinke  DNA/FCTT,  for  this  scoping  study  we  decided  to  represent  the  SHIST 
site  with  two  distinct  layers.  The  top  15  meters  were  modeled  with  a  seismic  velocity  of 
1 ,550  m/s  below  15  meters  was  a  layer  of  material  with  a  seismic  velocity  of  3,400  m/s. 
Since  there  is  no  experimental  data  to  indicate  the  stress  level  for  joint  “lock  up”,  simple 
models  were  developed  to  approximate  the  behavior  of  the  jointed  material. 
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Previous  work  on  jointed  material  [5]  has  indicated  that  the  joints  will  “lock  up” 
anywhere  from  0.5  to  1  kbar  after  which  the  material  behaves  as  though  jointing  did  not 
exist.  Since  we  don’t  have  any  experimental  data  for  jointing  characteristics  of  the 
SHIST  site,  it  was  assumed  that  the  joints  would  lock  up  at  approximately  0.5  kbar.  Two 
new  models  with  seismic  velocities  of  3,400  m/s  (material  model  GRAN5B,  which  in  this 
report  will  be  referred  to  as  SHIST2)  and  1,500  m/s  (material  model  GRAN4B,  which  in 
this  report  will  be  referred  to  as  SHIST1)  were  developed  using  the  material  model 
GENERIC3  as  a  starting  point.  The  initial  wave  speeds  were  lowered  to  the  appropriate 
values  and  the  porosity  increased  so  that  stress-strain  response  of  the  jointed  material 
models  would  intersect  with  the  baseline  model  at  0.5  kbar,  Figure  2-5. 


Figure  2-5.  Schematic  representation  of  elastic  stress-strain  assumption  used  in  the 
development  of  SHIST1  and  SHIST2  jointed  granite  material  models. 
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Since  it  is  assumed  that  the  material  is  elastic  for  stresses  below  0.5  kbar,  a 
closed  form  solution  for  this  intersection  is  possible.  The  constrained  modulus,  M,  can 
be  calculated  from  the  seismic  velocity,  c,  and  the  initial  density,  po,  by  the  relationship 

M  =  PoC^  (2.3) 

With  a  Poisson’s  ratio  of  0.25,  the  bulk  modulus,  K,  can  be  calculated  as 

_  M(1  +  "o)  _  M  _  PqC  /p 

Also,  for  uniaxial  strain,  the  stress,  a,  can  be  calculated  as  a  function  of  compression,  ji, 
or  density,  p,  by  the  relationship 

a  =  Mn  =  1.8K£^  (2.5) 

Po 

Using  the  model  GENERIC3  as  the  baseline  “NTS  like”  jointed  material  (seismic 
velocity  of  4,500  m/s  and  bulk  modulus  of  0.3  Mbar)  the  density  at  which  the  granite 
reaches  0.5  kbar  in  uniaxial  strain  Is  2.626  (or  0.09%  strain).  Assuming  that  the  SHIST 
jointed  granite  will  lock  up  at  this  same  stress  and  density,  that  it’s  stress  strain  behavior 
converges  with  model  GENERIC3,  the  initial  density  of  model  SHIST2  (seismic  velocity 
of  3,400  m/s  and  bulk  modulus  of  0.169  Mbar)  is  2.621  g/cc  and  the  Initial  density  of 
model  SHIST1  (seismic  velocity  of  1,550  m/s  and  bulk  modulus  of  0.035  Mbar)  is  2.605 
g/cc. 


The  SHIST  material  models  developed  by  this  elastic  approximation  are  intended 
only  for  use  in  the  current  SHIST  scoping  calculation  series.  The  in-situ  seismic 
velocities  used  here  were  obtained  from  low  stress  geophone  data,  and  there  are  large 
uncertainties  in  the  stress  level  at  which  the  joints  close.  If  the  SHIST  program  is 
reactivated  in  the  future,  it  is  recommended  that  additional  bulk  material  properties  data 
be  obtained  and  laboratory  tests  conducted  at  various  loading  rates  for  use  in  the 
development  of  more  representative  material  models. 
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The  properties  of  material  models  SHIST1,  SHIST2  and  GENERIC3  are 
summarized  in  Table  2-2  below.  For  reference  purposes,  the  SHEL  equation  of  sate 
input  parameters  have  been  included  in  Appendix  A  and  the  CRALE  input  decks  in 
Appendix  B. 


Table  2-2.  Properties  of  granite  material  models  SHIST1 ,  SHIST2,  and  GENERIC3. 


Model 

Density 

g/cc 

%AFV 

Initial 

K 

Mbar 

GENERIC1 

(GRAN1E) 

2.624 

2.65 

1.0 

0.3 

4,500 

SHIST2 

(GRAN5B) 

2.621 

2.65 

1.1 

0.169 

3,400 

>  15  m 

KSilUH 

2.605 

2.65 

1.7 

0.035 

1,550 

0-15m 

Uniaxial  stress-strain  response  curves  for  SHIST1  (GRAN4B)  and  SHIST2 
(GRAN5B)  are  shown  in  Figures  2-6a,b  and  2-7a,b.  In  each  of  these  figures, 

GENERIC3  (GRAN1 E)  is  shown  with  the  equivalent  strain  offset  used  in  the  calculations 
above  to  estimate  joint  lockup  at  0.5  kbar.  As  can  be  seen,  the  highly  jointed  SHIST 
materials  come  close  to  approximating  the  NTS  GENERIC3  granite  response  for 
stresses  above  0.5  kbar.  Since  these  “jointed”  models  were  intended  to  be  used  for  first 
order  scoping  calculations,  the  differences  observed  above  0.5  kbar  were  not 
considered  critical  at  this  time.  If  future  work  is  warranted,  they  will  be  modified  so  that 
a  better  match  with  the  generic  model  occurs  after  joint  lockup.  It  can  be  seen  in  Figure 
2-7  that  the  SHIST1  (GRAN4B)  model  does  not  converge  with  the  GENERICS  model 
until  stresses  exceed  1 .5  kbar.  This  is  due  to  the  differences  between  the  exponential 
nature  of  the  SHEL  EOS  and  the  linear  elastic  approximations  made  in  the  development 
of  GRAN4B.  These  differences  are  less  apparent  in  SHIST2  (GRAN5B,  1.1%  AFV) 
since  it  is  less  hysteretic  than  SHIST1  (GRAN4B,  1 .7%  AFV). 
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(a)  stress-strain  to  8  kbar 


Figure  2-6.  Uniaxial  stress  response  of  SHIST2  (GRAN5B) 
and  GENERIC3  (GRAN1E).  GENERIC3  is  also 
shown  with  a  0.1%  strain  offset. 
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(b)  stress-strain  to  80  kbar 


Figure  2-6.  Uniaxial  stress  response  of  SHIST2  (GRAN5B) 
and  GENERIC3  (GRAN1E).  GENERIC3  is  also 
shown  with  a  0.1%  strain  offset.  (Continued) 
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(a)  stress-strain  to  8  kbar 


Figure  2-7.  Uniaxial  stress  response  of  SHIST1  (GRAN4B)  and 
GENERIC3  (GRAN1 E).  GENERIC3  is  also  shown 
with  a  0.7%  strain  offset. 
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(b)  stress-strain  to  80  kbar 


Figure  2-7.  Uniaxial  stress  response  of  SHIST 1  (GRAN4B)  and 
GENERIC3(GRAN1  E).  GENERIC3  is  aiso  shown 
with  a  0.7%  strain  offset.  (Continued) 
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2.4  CONCRETE  STEMMING  MODEL  DESCRIPTION. 


The  stemming  material  used  in  this  study  is  HJC7  concrete  modeled  with 
the  MClST  material  model.  This  material  model  (CON5K4)  has  been  used  successfully 
in  the  past  (HUSKY  JAGUAR  and  MIGHTY  NORTH)  and  was  used  here  unchanged. 
The  uniaxial  stress-strain  response  of  CON5K4  is  compared  to  the  SHIST2  (GRAN5B) 
granite  model  in  Figure  2-8.  It  can  be  seen  that  the  concrete  model  is  similar  to  the 
SHIST  model  to  a  stress  level  of  approximately  1  kbar,  after  which  it  is  much  more 
hysteretic.  Therefore,  for  stresses  above  1  kbar,  the  shock  wave  in  the  granite  would 
be  expected  to  outrun  the  shock  in  the  concrete  filled  emplacement  hole.  This  effect 
should  be  localized  to  the  vicinity  of  the  emplacement  hole  and  would  not  be  expected 
to  Influence  the  free  field  environment  at  the  SHIST  experimental  gage  locations.  For 
reference,  the  CRALE  input  deck  for  the  CON5K4  concrete  model  has  been  included  in 
Appendix  B. 
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Figure  2-8.  Uniaxial  stress-strain  response  of  material  models  CON5K4 
concrete  and  SH1ST2  (GRAN5B)  jointed  granite. 
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SECTION  3 


PRELIMINARY  FREE-FIELD  ENVIRONMENT 
ID  CALCULATIONS  FOR  A  20-TON  CONTAINED  CHARGE 

(TASK  2) 


3.1  1 D  CALCULATION  SETUP. 

The  proposed  SHIST  experiment  was  designed  specifically  to  generate  a 
spherical  shock  wave  into  the  free-field  environment.  The  fully  contained  20-ton  QM- 
100/R  high  explosive  charge  was  to  be  placed  in  a  spherical  cavity  to  eliminate  any 
edge  effects  that  may  occur  from  a  cylindrical  charge  configuration.  The  experiment 
therefore  lends  itself  well  to  the  use  of  1 D  code  calculations  for  examination  of  the 
ground  shock  phenomenology. 

To  evaluate  the  shock  propagation  from  the  charge,  a  number  of  1-D  CRALE 
calculations  were  performed.  A  JWL  detonation  model  for  QM-1 00/R  was  used  in  all 
the  ID  (and  2D)  calculations.  This  model  was  developed  by  LLNL,  having  an  initial 
density  of  1 .509  g/cc  and  Chapman-Jouget  detonation  pressure  of  200  kbar.  For 
reference,  the  input  deck  for  the  QM-1 00/R  HE  EOS  used  in  the  CRALE  code  has  been 
included  in  Appendix  B. 

The  CRALE  1  -D  calculations  were  run  with  an  average  zone  size  of  20  cm  in  the 
granite  and  with  20  cells  in  the  HE  sphere  (radius  of  1 .42  m).  The  CRALE  HE  bum 
model  was  used  to  create  the  initial  shock  wave  into  the  surrounding  media.  Two  series 
of  calculations  were  run,  one  using  the  SHIST1  and  SHIST2  material  models  for  jointed 
granite  described  in  Section  2,  and  a  second  using  the  GENERIC2  and  GENERIC3 
granite  material  models  to  evaluate  potential  perturbations  due  to  material  model 
uncertainties. 
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3.2  1 D  CALCULATIONS  USING  SHIST1  AND  SHIST2  JOINTED  GRANITE 

MATERIAL  MODELS. 

Two  1D  calculations  were  run  with  the  SHIST1  (GRAN4B,  1,550  m/s, 
1.7%  AFV)  and  SHIST2  (GRAN5B,  3,400  m/s,  1.1%  AFV)  “jointed”  granite  material 
models  described  in  Section  2.3  above.  These  calculations  were  run  with  an  average 
zone  size  of  20  cm  in  the  granite  and  with  20  cells  in  the  QM-100/R  HE  sphere  (radius 
of  1 .42  m).  Figures  3-1  through  3-3  show  the  peak  stress,  velocity,  and  displacement 
vs.  range  attenuation  curves  calculated  for  the  two  SHIST  material  models.  The 
calculated  stress,  velocity  and  displacement  waveforms  for  SHIST1  and  SHIST2  are 
plotted  at  ranges  between  2  and  100  meters  in  Figures  3-4  through  3-6. 
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Figure  3-1 .  Calculated  1 D  peak  stress  vs.  range  for  a  20-ton  QM-1 00/R 
charge  in  SHIST1  (GRAN4B,  1 ,550  m/s.  1.7%  AFV  and 
SHIST2  (GRAN5B,  3,400  m/s,  1.1%  AFV)  jointed  granite. 
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Figure  3-2  Calculated  1 D  peak  velocity  vs.  range  for  a  20-ton  QM-100/R 

charge  in  SHIST1  (GRAN4B.  1 ,550  m/s,  1.7%  AFV  and  SH1ST2 
(GRAN5B,  3,400  m/s,  1.1%  AFV)  jointed  granite. 
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Figure  3-3.  Calculated  1 D  peak  displacement  vs.  range  for  a 
20-ton  QM-100/R  charge  in  SHIST1  (GRAN4B, 
1,550  m/s,  1.7%  AFV  and  SHIST2  (GRAN5B, 
3,400  m/s,  1.1%  AFV)  jointed  granite. 
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Figure  3-4.  Calculated  1 D  stress  waveforms  for  a  20-ton  QM-1 00/R 
charge  in  SHIST1  (GRAN4B.  1,550  m/s,  1.7%  AFV  and 
SHIST2  (GRAN5B,  3,400  m/s,  1.1%  AFV)  jointed  granite 
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Figure  3-5.  Calculated  1 D  velocity  waveforms  for  a  20-ton  QM-1 00/R 
charge  in  SHIST1  (GRAN4B,  1,550  m/s,  1.7%  AFV  and 
SHIST2  {GRAN5B,  3,400  m/s,  1.1%  AFV)  jointed  granite. 
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Figure  3-6.  Calculated  ID  displacement  waveforms  for  a  20-ton  QM-100/R 
charge  in  SHIST1  (GRAN4B,  1,550  m/s,  1.7%  AFV  and  SHIST2 
(GRAN5B,  3,400  m/s,  1.1%  AFV)  jointed  granite. 
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As  can  be  seen  in  the  figures  above,  the  1-D  calculation  for  material  model 
SHIST2  (GRAN5B,  3,400  m/s,  1 .1%  AFV)  shows  a  distinct  range  at  which  the  onset  of 
the  “elastic”  attenuation  regime  begins,  approximately  30  meters  from  the  charge  center 
at  a  stress  level  of  0.3  kbar.  This  onset  occurs  at  the  point  where  the  peak  stress  and 
velocity  attenuation  rates  approach  a  slope  of  1  (the  theoretical  spherical  attenuation  of 
a  purely  elastic  wave).  As  can  be  seen  in  Figures  3-1  and  3-2  above,  the  SHIST2 
plastic  wave  has  a  plastic  stress  wave  attenuation  on  the  order  of  a  plastic  velocity 
wave  attenuation  on  the  order  of  R'^  ®,  and  an  elastic  wave  attenuation  of  R  .  Since 
most  of  the  experimental  gages  will  be  placed  within  the  SHIST2  layer  (depth  >  15  m), 
these  1 D  attenuations  can  be  used  to  assist  In  gage  placement.  Potential  surface 
rarefaction  will  be  examined  in  the  2-D  calculations  reported  in  Section  4  below. 


The  1-D  results  using  material  model  SHIST1  (GRAN4B,  1,550  m/s,  1.7%  AFV) 
does  not  show  a  distinct  range  at  which  the  shock  attenuates  as  an  elastic  wave. 
SHIST2  is  has  a  lower  seismic  velocity  and  is  more  hysteretic  than  SHIST1  which 
causes  the  shock  to  attenuate  more  rapidly.  An  “elastic”  regime  is  not  predicted  by  the 
calculation  within  a  range  of  100  meters.  Figures  3-1  and  3-2.  Since  this  material  model 
is  a  oversimplification  of  the  weathered  rock,  it  may  be  an  overestimate  of  the  jointing 
effects  of  the  upper  layer  but  without  experimental  data,  this  is  difficult  to  confirm. 

3.3  ID  CALCULATIONS  USING  GENERIC2  AND  GENERIC3  GRANITE 

MATERIAL  MODELS. 

Figures  3-7  through  3-9  show  the  peak  stress,  velocity,  and  displacement  vs. 
range  attenuation  curves  calculated  for  the  GENERIC2  (GRAN1C,  4,500  m/s,  1%  AFV) 
and  GENERIC3  (GRAN1E,  4,500  m/s,  1%  AFV)  granite  material  models.  These 
calculations  were  included  to  examine  material  properties  excursions  to  the  SHIST  site. 
Since  accurate  in-situ  data  is  not  currently  available  for  the  SHIST  site,  the  degree  of 
weathering  in  uncertain  and  it  Is  not  known  at  what  stress  level  the  joints  will  close.  The 
GENERIC2,3  models  are  a  higher  quality  granite  than  the  SHIST1,2  models  and  was 
felt  their  predicted  attenuations  could  assist  in  future  test  gage  planning.  As  can  be 
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seen,  the  onset  of  the  “elastic”  attenuation  regime  for  the  GENERIC  models  occurs  at  a 
slightly  greater  range  than  predicted  for  SHIST1,  approximately  40  meters  from  the 
charge  center  at  a  stress  level  of  0.3  kbar.  The  calculated  stress,  velocity  and 
dispiacement  waveforms  for  GENERiC2  and  GENERiCS  are  plotted  at  ranges  between 
2  and  100  meters  in  Figures  3-10  through  3-12.  The  GENERIC2,3  models  predict  a 
plastic  stress  wave  attenuation  on  the  order  of  R’’  ®,  a  piastic  velocity  wave  attenuation 
on  the  order  of  R’’  and  an  elastic  wave  attenuation  of  R*^  \ 

The  two  GENERIC  1%  AFV  granite  modeis  can  be  seen  to  give  virtualiy  identicai 
shock  attenuations.  Slight  differences  can  be  seen  due  to  variations  in  the  strength 
characteristics  of  the  two  materials  (linear  DPA/M  yieid  surface  vs.  exponential  surface). 
GENERIC3  gives  slightly  smaller  displacements  than  GENERIC2  but  times  of  arrival  are 
identical  in  both  materials,  with  the  elastic  wave  running  at  4,500  m/s. 


33 


RANGE  (M) 


Figure  3-7.  Calculated  1 D  peak  stress  vs.  range  for  a  20-ton  QM-100/R 
charge  in  4,500  m/s,  1%  AFV  granite  using  GENERIC2 
(GRAN1C)  and  GENERIC3  (GRAN1E)  material  models. 
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Figure  3-8.  Calculated  1 D  peak  velocity  vs.  range  for  a  20-ton  QM-1 00/R 
charge  in  4,500  m/s,  1%  AFV  granite  using  GENERIC2 
(GRAN1C)  and  GENERIC3  (GRAN1E)  material  models. 


Figure  3-9.  Calculated  1 D  peak  displacement  vs.  range  for  a  20-ton  QM-100/R 
charge  in  4,500  m/s,  1%  AFV  granite  using  GENERIC2  (GRAN1C) 
and  GENERIC3  (GRAN1E)  material  models. 
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Figure  3-10.  Calculated  ID  stress  waveforms  range  for  a  20-ton  QM-100/R 

charge  in  4,500  m/s,  1%  AFV  granite  using  GENERIC2  (GRAN  1C) 
and  GENERIC3  (GRAN1 E)  material  models. 
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Figure  3-12.  Calculated  ID  displacement  waveforms  range  for  a  20-ton  QM-100/R 
charge  in  4,500  m/s,  1%  AFV  granite  using  GENERIC2  (GRAN  1C)  and 
GENERIC3  (GRAN1E)  material  models. 


SECTION  4 


2D  SCOPING  CALCULATIONS  FOR  SHIST 
(TASK  3) 


4.1  2D  CALCULATION  SETUP. 

The  material  models  discussed  in  the  sections  above  were  implemented  in  the 
CRALE  2-D  code  and  used  to  investigate  the  details  of  the  interaction  of  the  explosively- 
driven  shock  with  the  surrounding  concrete  stemming  materials  and  granite.  We  used 
the  latest  proposed  charge  and  stemming  configuration  provided  to  us  by  DNAFC.  This 
plan  called  for  a  20  ton  sphere  of  QM-100/R  with  a  diameter  of  2.896  meter  (9.5  ft)  to  be 
placed  at  a  shot  depth  of  35  meters.  A  1 .066  meter  (42  inch)  diameter  downhole  stem 
extends  vertically  from  the  surface  to  approximately  15  meters  below  charge  (50  meter 
depth).  The  downhole  was  to  be  filled  with  HJC7  concrete  for  which  an  existing  CRALE 
material  model  exists  (CON5K4  using  the  MClST  EOS). 

Three  2D  calculations  were  performed  in  this  study  to  examine  the  shock 
propagation  in  the  granite  materials  discussed  in  the  sections  above,  and  to  evaluate 
potential  perturbations  due  to  the  material  model  uncertainties  of  the  SHIST  site.  Each 
calculation  was  run  to  a  maximum  time  of  20  msec  to  capture  the  peak  displacements  at 
the  experimental  gage  locations.  The  first  calculation,  designated  2DSHIST03  and 
discussed  In  Section  4.2,  consisted  of  a  top  layer  of  SHIST1  (GRAN4B,  1,550  m/s,  1.7% 
AFV)  jointed  granite  to  a  depth  of  15  meters  below  which  was  an  infinite  layer  of  SHIST2 
(GRAN5B,  3,400  m/s,  1 .1%  AFV)  jointed  granite.  This  calculation  Is  currently  our  best 
representation  for  the  SHIST  experiment.  The  second  calculation,  designated 
2DSHIST02  and  discussed  in  Section  4.3,  consisted  of  a  homogeneous  layer  of  SHIST2 
jointed  granite.  This  calculation  could  be  used  in  conjunction  with  2DSHIST03  to  help 
evaluate  the  potential  layering  effects  on  the  free  field  waveforms.  The  third  calculation, 
designated  2DSHIST01  and  discussed  in  Section  4.4,  consisted  of  a  top  layer  of 
GENERIC2  (GRAN1C,  4,500  m/s,  1%  AFV)  granite  below  which  was  an  Infinite  layer  of 
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GENERICS  (GRAN1E,  4,500  m/s,  1%  AFV)  granite.  This  calculation  was  used  as  a 
second  material  model  excursion  and  2D  zoning  study  to  be  compared  to  the  1D 
calculations  discussed  in  Section  3. 

Due  to  zoning  constraints  of  the  2D  code,  the  calculations  were  started  with  an 
initial  zone  size  of  approximately  10  cm  in  the  vicinity  of  the  HE  charge.  The  initial  2D 
calculation  grid  is  shown  In  Figure  4-1.  The  HE  was  modeled  as  a  spherical  charge  with 
20  cells  in  the  radial  direction  and  the  CRALE  2D  bum  logic  was  used  to  calculate  its 
detonation.  After  the  HE  was  complete  and  the  shock  wave  entered  the  surrounding 
media,  at  a  time  0.4  msec,  the  grid  was  dezoned  to  an  average  zone  size  of  40  cm. 

This  zone  size  was  then  maintained  throughout  the  remainder  of  the  calculation,  roughly 
twice  that  used  in  the  1 D  calculations  reported  in  Section  3. 

An  array  of  1 97  target  points  were  placed  in  the  2D  calculations  to  a  range  and 
depth  of  100  meters.  Included  in  this  array  were  target  points  placed  at  the  latest 
experimental  gage  locations  that  you  provided  us.  Figure  4-2  shows  the  target  point 
array,  with  the  experimental  gages  bolded. 
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Figure  4-1 .  Initial  grid  used  in  2D  CRALE  SHIST  scoping  calculations. 


15. 
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Figure  4-2.  Target  point  locations  and  site  layering  for  2D  CRALE 
SHIST  scoping  calculations. 
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4.2  CRALE  2D  CALCULATION  2DSHiST03. 

SHIST1  (GRAN4B,  1,550  m/s,  1.1%  AFV)  0  to  15  meter 

SHiST2  (GRAN5B,  3,400  m/s)  depth  >  15  meters 

The  first  calculation  of  this  study,  2DSHIST03,  was  run  on  the  DNA  Field 
Command  Silicon  Graphics  Crimson  workstation.  CRALE  has  been  successfully  ported 
to  this  platform  and  the  code  is  fully  operational  there  (the  current  version  of  CRALE  is 
CRALEPL3).  The  2DSHIST03  calculation  examined  the  2D  shock  wave  propagation 
through  a  layered,  highly  weathered  and  jointed  granite  site.  The  calculation  consisted 
of  material  model  SHIST1  (GRAN4B,  1 ,550  m/s,  1 .7%  AFV)  to  a  depth  of  15  meters 
with  model  SHIST2  (GRAN5B,  3,400  m/s,  1.1%  AFV)  for  depths  greater  than  15  meters. 
In  this  calculation,  the  shock  in  layer  SHIST2  can  be  considered  spherical  (ID)  up  to  the 
time  it  hits  the  15  meter  deep  layer  interface.  The  presence  of  the  concrete  filled 
emplacement  hole  does  not  appear  to  affect  the  behavior  of  the  free-field  environment 
in  the  region  of  the  experimental  gages.  Once  the  shock  enters  the  SHIST1  layer,  it 
begins  to  propagate  at  a  slower  velocity  through  the  upper  layer.  This  results  in  the 
formation  of  a  downward  directed  relief  wave  into  the  SHIST2  layer.  A  second  surface 
rarefaction  wave  results  when  the  shock  reaches  the  ground  surface.  Figures  4-3 
through  4-5  show  velocity  vector  fields  from  the  calculation  at  8, 16,  and  22  msec 
respectively  in  which  the  surface  relief  as  well  as  relief  from  the  15  meter  interface  can 
be  observed.  The  vertical  and  horizontal  stress,  velocity,  and  displacement  waveforms 
for  the  target  points  at  the  experimental  gage  locations  are  included  as  Figures  4-6 
through  4-25.  The  peak  stress  contours  for  the  2DSHIST03  calculation  are  shown  in 
Figure  4-26. 
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Figure  4-3.  Velocity  vector  field  at  8  msec  for  2DSHIST03  CRALE  calculation. 
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Rgure  4-4.  Velocity  vector  field  at  1 6  msec  for  2DSHIST03  CRALE  calculation. 
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Figure  4-5.  Velocity  vector  field  at  22  msec  for  2DSHIST03  CRALE  calculation. 
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Figure  4-12.  Target  point  35  time  history  for  2DSHIST03  CRALE  calculation. 
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Figure  4-1 5.  Target  point  44  time  hi 
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Figure  4-1 6.  Target  point  46  time  history  for  2DSHIST03  CRALE  calculation. 
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Figure  4-20.  Target  point  70  time  history  for  2DSHIST03  CRALE  calculation 


62 


VELOCITY(M/S)  DISPLflCE(M)  STRESS+Q( KBHR ) 

crZ»00 _ 0; _ 2.00  -Q.Ql _ 0^ _ Q.Ql  0- _ 0*20 _ 0 


SHIST03  20TON  QHIOOR  IN  0RRN6B  BNO  0RnN4B 


17IZ9IS3 


17-Rtt9-9 

16. 


95. 

RfWOC  CH)  «  20.0 
DEPTH  IH)  s-50.0 


X  s  HaRIZeNTRL 
Y  s  VERTICAL 


Q.Q05  0.010  0.015  0.020  0.025  0.030 

TIMECSEC) 


Figure  4-23.  Target  point  95  time  history  for  2DSHIST03  CRALE  calculation. 


Figure  4-24.  Target  point  96  time  history  for  2DSHIST03  CRALE  calculation. 
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SHIST03  20  TON  QM100R  IN  GRAN5B/GRAN4B 
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Figure  4-26.  Peak  stress  contours  for  the  2DSHIST03  CRALE  calculation 
SHIST1  (GRAN4B,  1,550m/s,  1.7%AFV)  to  15  meter  depth 
with  SHIST2  (GRAN5B,  3,400m/s,  1.1%AFV)  below. 
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4.3  CRALE  2D  CALCULATION  2DSHIST02. 


HOMOGENEOUS  SHIST2  (GRAN5B,  3,400  M/S,  1.1%  AFV) 

The  second  calculation  of  this  study,  2DSHIST02,  consisted  of  a  homogeneous 
site  of  SHIST2  (GRAN5B,  3,400  m/s,  1.1%  AFV).  The  shock  propagation  in  this 
calculation  is  spherical  until  it  hits  the  free  surface,  at  which  time  a  relief  wave 
propagates  back  into  the  ground.  Figures  4-27  through  4-29  show  velocity  vector  fields 
from  the  calculation  at  8, 1 6,  and  22  msec  respectively  in  which  the  surface  relief  can  be 
observed.  The  vertical  and  horizontal  stress,  velocity,  and  displacement  waveforms  for 
the  target  points  at  the  experimental  gage  locations  are  included  as  Figures  4-30 
through  4-50.  The  peak  stress  contours  for  the  2DSHIST02  calculation  are  shown  in 
Figure  4-51 . 
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Figure  4-27.  Velocity  vector  field  at  8  msec  for  2DSHIST02  CRALE  calculation. 
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Figure  4-28.  Velocity  vector  field  at  1 6  msec  for  2DSHIST02  CRALE  calculation 
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Figure  4-29.  Velocity  vector  field  at  22  msec  for  2DSHIST02  CRALE  calculation. 
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SHIST02  20TeN  OtllOOR  IN  GRflNSB 


Figure  4-31 .  Target  point  4  time  history  for  2DSHIST02  CRALE  calculation 
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Figure  4-36.  Target  point  33  time  history  for  2DSHIST02  CRALE  calculation. 
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SHIST02  20T0N  OHIOOR  IN  ORRNSB 


Figure  4-37.  Target  point  35  time  history  for  2DSHIST02  CRALE  calculation. 
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SHIST02  20TeN  OHIOOR  IN  GRRN5B 


Figure  4-39.  Target  point  43  time  history  for  2DSHIST02  CRALE  calculation. 
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SHIST02  20TON  OMIOOR  IN  ORf)N5B 


Figure  4-40.  Target  point  44  time  history  for  2DSHIST02  CRALE  calculation. 
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8HI8T02  20T8N  QniOOR  IN  ORRN5B 


Figure  4-42.  Target  point  48  time  history  for  2DSH1ST02  CRALE  calculation. 


SHIST02  20TeN  OHIOOR  IN  ORNNSB 


Figure  4-43.  Target  point  68  time  history  for  2DSHIST02  CRALE  calculation  . 


86 


8HIST02  20TQN  QHIOOR  IN  ORRNSB 


Figure  4-44.  Target  point  69  time  history  for  2DSHIST02  CRALE  calculation. 
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SH1ST02  ZOTON  QHIOOR  IN  GRANSB 


Figure  4-45.  Target  point  70  time  history  for  2DSHIST02  CRALE  calculation. 
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Figure  4-46.  Target  point  72  time  history  for  2DSHIST02  CRALE  calculation. 
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8HIST02  20TeN  OHIOOR  IN  ORANSB 


Figure  4-47.  Target  point  74  time  history  for  2DSHIST02  CRALE  calculation. 
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SHIST02  20TON  QHIOOR  IN  ORRNSB 


Figure  4-48.  Target  point  95  time  history  for  2DSH1ST02  CRALE  calculation. 
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Figure  4-50.  Target  point  108  time  history  for  2DSHIST02  CRALE  calculation. 
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Figure  4-51 .  Peak  stress  contours  for  the  2DSHIST02  CRALE  calculation 
Homogeneous  SHIST2  (GRAN5B,  3,400mps,  1.1%AFV). 
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4.4  CRALE  2D  CALCULATION  2DSHIST01. 

GENERIC2  (GRAN1 C,  4,500  m/s)  0  to  1 5  meter 

GENERIC3  (GRAN1 E,  4500  m/s)  depth  >  15  meters 

The  2DSH1ST01  calculation  was  run  with  material  model  GENERIC2  (GRAN1C) 
in  the  top  15  meters  and  material  model  GENERIC3  (GRAN1E)  in  the  bottom  layer. 
Since  the  seismic  velocities  of  GENERIC1  and  GENERIC2  (4,500  m/s)  are  higher  than 
those  expected  in  the  field,  this  calculation  was  run  primarily  as  a  zoning  study  to 
compare  with  the  1 D  calculations  run  earlier.  Also,  since  GENERIC2,3  are  virtually 
identical  materials,  the  site  profile  of  this  calculation  can  be  considered  homogeneous 
and  the  shock  propagation  is  spherical  until  it  hits  the  free  surface.  The  2D  velocity  and 
displacement  waveforms  from  target  points  along  the  35  meter  depth  radial  are 
compared  to  the  GENERIC3  (GRAN1E)  ID  calculation  in  Figures  4-52  and  4-53.  As 
can  be  seen,  since  the  zoning  in  the  2D  calculation  was  coarser  than  the  1 D,  the  rise 
time  of  the  waveforms  are  longer.  However,  it  can  also  be  seen  in  the  displacement 
waveforms,  that  impulse  has  been  conserved  and  the  displacements  predicted  in  the 
SHST01  2D  calculation  are  equivalent  to  those  in  the  1 D  calculation.  The  peak  stress 
contours  from  the  SHIST01  2D  calculation  are  shown  in  Figure  4-54. 
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Figure  4-53.  35  meter  depth  radial  displacements  from  2DSHIST01 
calculation  compared  to  the  GENERIC3  (GRAN1E) 

1D  calculation. 
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SHIST01  20  TON  QM100R  IN  GRAN1E/GRAN1C 
TIME=  .020644  SEC  CYCLE  =  1752 


Figure  4-54.  Peak  stress  contours  for  the  2DSHIST01  calculation. 

GENERIC2  (GRAN  1C,  4.500m/s,  1%AFV)  to  15  meter 
depth  with  GENERIC3  (GRAN1E,  4500m/s,  1%AFV)  below. 
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APPENDIX  A 


SHELL  EQUATION  OF  STATE  PARAMETERS 


A.1  GENERAL 

Given:  a  state  of  the  material  defined  by  its  current 

a)  Density,  p(g/cc)  and  strain  state,  Ejj; 

b)  Specific  internal  energy,  E  (Terg/g) 

Find:  the  pressure  or  mean  stress,  P  (Mb)  and  the  deviatoric  stress  tensor,  a'jj 
Definitions: 

pQ,  the  initial  density 

P|.^  the  reference  density,  i.e.,  the  grain  density  of  the  solid 
Compression,  11  =  p/p^ 

(Note:  compression  is  relative  to  the  initial  density) 

Excess  compression,  |i  =  n  - 1 
Porosity  (total  void  volume),  v^  =  1  -  p^/p^ 
volumetric  strain,  e  =  -p/Ti 

Then: 


P(p,E)  =  P,(p,E)  +  P.(p.E)  (A.1) 

P,  (A.2) 

P.  =  rpE';  (E>E„)  (A.3) 

where 

|i'=p  +  PE  (A.4) 

E'  =  (E-E„)(1-e-'^-^l'^)  (A.5) 


A-1 


1  r  =  Y  - 1  =  0.6  +  0.023[ln(p  /  E)f  +  0.1 4  •  ln(p  /  E)  +  0.05  •  In  p  (A.6) 

Five  parameters  are  required  in  Equations  A.2  through  A.6  to  compute  the  pressure  in  a 
material  with  no  porosity  (voids),  i.e.  pQ=pf,  in  addition  to  the  reference/initial  density. 
They  are: 

=  minimum  energy  required  to  activate  P^,  (Terg/g) 

Bq  =  the  zero  pressure  bulk  modulus,  dP/dp  (Mb) 

=  the  maximum  bulk  modulus,  (Mb) 
p*  =  the  excess  compression  at  which  the  bulk  modulus  has 
increased  exponentially  to  ~[2/3  B^+1/3  B^] 

P  =coefficientof  thermal  expansion,  (g/Terg); 

BmPE  is  equivalent  to  the  Gruneisen  gamma  term  (TpE) 
for  the  energy  dependence  of  the  solid  pressure. 

A.2  HYSTERESIS  DUE  TO  AIR-FILLED  POROSITY. 

For  a  material  which  contains  air-filled  porosity,  i.e.,  less  than  100%  saturation, 
the  pores  may  collapse  irreversibly  under  compression  so  that  the  behavior  of  the 
material  is  path  dependent.  Several  additional  steps  and  a  new  independent  variable 
are  required  to  compute  the  pressure.  In  the  SHEL  routine,  the  new  variable,  pmax 
(Pj^),  is  the  maximum  excess  compression  ever  attained  by  the  material.  For  a  material 
containing  air-filled  void  spaces,  the  initial  density,  p^,  is  always  less  then  the  reference 
density,  i.e.;  Pq<p^..  By  convention,  the  excess  compression,  p,  relative  to  p^  is 
passed  to  the  SHEL  routine.  Since  the  SHEL  EOS  algorithms  are  based  on  the 
reference  density,  p^.,  the  routine  transforms  the  p's  to  correspond  to  p^,. 


^The  multiplication  symbol,  •,  has  been  included  where  needed  for  clarity:  In  is  the  natural  log;  and  7  is 
the  gamma-law  gas  EOS  parameter  (not  explicitly  used). 
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A  fundamental  assumption  governs  the  pressure  behavior  of  the  porous  material 
under  loading,  namely,  there  is  a  transformation  in  \i  to  an  equivalent  value  along  the 
load  path  of  the  solid  (void-free)  material.  To  calculate  the  pressure,  we  first  find  the 
solid  |i  equivalent  to  and  then  obtain  a  modified  p'  for  the  current  compressional 
state  of  the  material  to  use  in  Equation  A.2  as  follows. 

Define: 


Tlz=Po/Pr.  Pz='nz-1^0.  ez=-Pz/Tlz  (A-7) 

Any  excess  compression,  p  relative  to  can  be  transformed  into  an  equivalent  p^. 
relative  to  p^.,  by 

Pr=TlzP  +  Pz  (A.8) 


As  for  a  nonporous  solid  (Equation  B.4),  we  again  augment  p  and  p^^  by  adding  the 
thermal  contribution  to  both; 


p"=p+PE 

Pm=Pm+PE 


(A.9) 


An  offset  parameter.  A,  is  computed  next  to  transform  the  compression  into  the  frame  of 
the  reference  density; 


A  =  P2(1  -  )  +  ^,(1  _  )e“2^im/ez  (A.1 0) 

where  a.| ,  a2  and  p.^  are  additional  input  parameters.  The  first  term  on  the  right-hand 
side  of  Equation  B.10  provides  a  smooth  convex  curve  starting  at  p=Pq  which 
asymptotically  approaches  the  reference  P-p  hydrostat;  the  second  term  adds  a  toe  to 
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the  initial  loading.  Physically,  the  first  term  represents  the  irreversible  pore  collapse  and 
the  breakdown  of  cementation  bonds  in  the  matrix,  while  the  second  term  is  a  measure 
of  the  initial  elastic  response  of  those  bonds. 

Pg  is  now  computed  as  in  Equation  B.2: 


where 


(A.11) 


p  =  ri2ti"+A  (A.12) 

and 

p*'=  |i*/[Tiz(1-a2e“^^'^^^^)  +  aie"“''^"’^^^]  (A.13) 

The  assumption  that  unloading  paths  for  partially  collapsed  porous  media  parallel 
the  solid  load-unload  path  has  been  used  for  many  earth  media,  particularly  when  load- 
unload  data  were  not  available  or  inconsistent.  Recommended  load-unload  paths  for 
the  porous  shales  provided  by  WES  did  not  conform  to  this  behavior.  Unloading  from 
peak  pressures  less  than  the  elastic  toe  were  elastic;  as  the  peak  pressure  increased 
above  the  crush  level  so  did  the  initial  unload  moduli.  Thus  the  unloading  paths  appear 
to  stiffen  and  fan  out  as  the  peak  pressure  increases.  This  effect  was  modeled  by 
simply  adjusting  the  effective  maximum  excess  compression  as  the  material  unloaded, 

i.e.,  forpxi<^m  <V-X2 


(A.14) 


where 


^  “  (M-m  M’X'|)^(M’X2  M^xi) 


(A.15) 


and  \Lx^and\ix2  are  additional  input  parameters. 

Thus,  for  a  porous  solid,  six  parameters  must  be  specified,  namely; 
pQ  =  the  initial  density  of  the  material  (Tipj.) 

a.j  =  relative  slope  of  initial  loading  on  the  elastic  toe  to  that 
of  the  non-porous  solid 

|i1  =  amount  of  compression  on  loading  toe 
02  =  a  fitting  parameter  that  governs  the  elastic  slope  of  the 
load  curve  while  crushing  out  the  voids 
p.x1  =  minimum  compression  to  begin  fanning  the  unload  paths 
ixx2  =  compression  at  which  unload  path  parallels  zero  porosity 
load-unload  curve. 
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A.3  SOLID-SOLID  PHASE  CHANGE. 


Many  materials  exhibit  a  solid-solid  phase  change  under  compression.  Based  on 
Hugoniot  and  release  data  for  silicates,  it  appears  that  under  shock  loading  such  a 
change  is  hysteretic,  i.e.,  on  unloading  the  material  initially  remains  in  the  higher  density 
state  and  does  not  revert  to  the  original  state  until  well  down  the  release  path.  Such 
behavior  can  be  modeled  by  assuming  that  the  phase  change  is  a  function  of  both 
compression  and  energy. 

To  compute  the  new  solid  P  ,  Equation  A.2  or  A.1 1  is  again  modified  by 

o 

replacing  the  constant  and  the  excess  compression,  p'  (Equation  A.4)  or  p 
(Equation  A.1 2),  by 


B'n,  and  p'-A,  or  p- Ap  (A.1 6) 

where 


Ap  = 


0 

g^2 


pE<(pE)^-  (all  state  1) 
pE>(pE), 


(A.17) 


I  Bm  pEs(pE)^:  (all  state  1)  (A.1 8) 

I  Bm+g(Bm2-Bm)  pE>(pE)i 


P2=p2/Pr-1. 


1  _  e-a3lpB-(pB)ll 

■j  +  g-a3[pE-(pE)il 


(A.19) 
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and  Bm2,  (pE)i ,  P2  and  as  are  material  constants.  Thus,  four  additional  input 
parameters  are  required  for  the  phase  change: 

Pg  =  the  reference  density  of  the  2nd  phase 
(pE).|  =  the  energy  per  unit  volume  at  which  the  phase 
change  begins 

Og  =  a  parameter  setting  the  rate  of  conversion  between 
the  phases 

®m2  “  maximum  bulk  modulus  (dp/dp)  of  the  2nd  phase. 

NOTE:  B^2  relative  to  p^,  of  the  1st  phase; 
the  modulus  relative  to  P2  is  =  (P2/Pf)  B^2 

A.4  DEVIATORIC  BEHAVIOR. 

The  complete  description  of  a  material's  stress-strain  behavior  requires  the 
calculation  of  the  deviatoric  stress  tensor, a'y,  in  addition  to  the  mean  stress  or  pressure. 

The  SHEL  EOS  model  assumes  that  the  stress  deviators  depend  on  the  strain 
deviators,  e'jj,  and  a  yield  function,  Y.  The  stresses  are  calculated  in  two  steps.  First, 

the  change  in  the  deviators  is  assumed  to  be  elastic  and  an  incremental  change  is 
added  to  the  previous  deviators, 

a'jj  (n)  =  a’jj  (n  - 1)  +  2Gde'|j  (A.20) 

where  G  is  the  shear  modulus  of  the  material  and  de'ij  is  the  change  in  the  ij-th 

component  of  the  deviatoric  strain  tensor  between  times  n-1  and  n.  Furthermore,  it  is 
assumed  the  shear  modulus  can  be  derived  from  the  current  bulk  modulus  through 
Poisson's  ratio  and  the  standard  elastic  relation 

2G  =  3B(1-2v)/(1  +  v)  (A.21) 


A-7 


After  the  stresses  have  been  updated  by  Equation  A.20  a  check  is  made  to  see  if  the 
elastic  changes  are  valid.  The  second  invariant  of  the  deviatoric  stress  tensor,  J'2  is 
compared  to  the  square  of  a  yield  function,  Y  which  depends  only  on  the  mean  stress, 
P.  If  J'2  exceeds  y2,  the  material  is  deforming  plastically  and  the  deviators  are 
proportionally  reduced  so  that  J'2  is  equal  to  y2. 

For  the  rock  media  of  interest,  two  yield  surface  functions  are  available,  an  exponential 
and  a  simple  Drucker-Prager  /  Von-Mises  limit  with  the  input  variable  Py  used  as  a  flag. 
When  1/Py  >10,  the  yield  function  is 


Y  =  max 


—P/P 


'  vm 


(A.22) 


The  four  input  parameters  required  to  calculate  the  deviatoric  stress  tensor  are; 
v,  Poisson's  ratio 

Yvm.  the  limiting  yield  strength  (kb) 

Yq,  the  yield  strength  at  zero  pressure  (kb) 

P  the  current  pressure 

Py  the  pressure  at  which  the  yield  function  -(2/3  Ym  + 1/3  Yq) 

When  1/Py  <10. 


where 


Y  =  max-^ 


Py 


Y 


vm 


Yvm,  the  limiting  yield  strength 

Yq,  the  yield  strength  at  zero  pressure  (cohesion) 

P  the  current  pressure 

1/Py  the  Drucker-Prager  slope 


(A.24) 
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Table  A-1.  SHELEOS  input  parameters. 


Value 


GENERIC2  I  GENERIC3  SHIST1 


GRAN  IE 


2.624 


GRAN4B 


2.605 


SHIST2 


GRAN5B 


2.621 


Comments 


Initial  density 


Reference  density 


Min  energy  for  Py 


Tensile  Limit 


initial  bulk  modulus 


Maximum  bulk 
modulus 


Rate  of  change  of  B 


End  of  elastic  toe 


Rate  of  increase  of  B 
on  elastic  toe 


.3 

i 

.3 

I 

Mbar  0.0125 


Mbar  0.0001 


Mbar  0.30 


0.007 


2.0 


18.7 


1.0 


3.0 


0.25 


0.0125 


0.0001 


30.0 


0.007 


2.0 


2.3 


1.7 


3.0 


0.25 


0.0125 


0.0001 


30.0 


0.0 


0.02 


0.09 


.3 


0.007 


2.0 


2.3 


1.1 


3.0 


0.25 


0.0125 


0.0001 


30.0 


Min  EN  for  solid 
in  2nd  phase 


Max  B  of 
2nd  phase  solid 


Porosity  (total 
void  volume) 


Thermal 
expansion  coef. 


Poisson’s  ratio 


Von  Mises  yield  limit 


Cohesion 


Slope  Coefficient 

<10  linear,  >10  exponential 
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APPENDIX  B 

CRALE  EOS  INPUT  DECKS 


This  Appendix  included  the  CRALE  input  decks  necessary  to  run  the  calculations 
discussed  in  this  report.  The  CRALE  EOS  input  parameters  are  discussed  further  in 
Reference  3. 


B.1  SHEL  EOS  MODELS  SHIST1  (GRAN4B  AT  1 ,550  M/S)  AND  SHIST2 
(GRAN5B  AT  3,400  M/S). 


123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789- 


GRAN4B  1  1 

1.7%  GRAN 

1550  m/sec 

EXP  yield  surface 

2.605 

2.650 

0.04 

0.0 

-120.E-06 

0.30 

0.7 

.03 

0.0 

0.117 

0.0 

1.00 

0.0 

0.020 

0.09 

0.0 

0.0 

0.0 

0.0 

0.05 

1.0 

3.0 

0.00 

0.0 

4.30 

.007 

.13 

2.0 

2.3 

0.25 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0125 

0.0001 

30.0 

0.0 

0.0 

0.0 

0.0 

0.0 

GRAN5B  1  1 

1.1%  GRAN 

3400  m/sec 

EXP  yield  surface 

2.621 

2.650 

0.04 

0.0 

-120.E-06 

0.30 

0.7 

.03 

0.0 

0.561 

0.0 

1.00 

0.0 

0.020 

0.09 

0.0 

0.0 

0.0 

0.0 

0.05 

1.0 

3.0 

0.00 

0.0 

4.30 

.007 

.13 

2.0 

2.3 

0.25 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0,0125 

0.0001 

30.0 

0.0 

0.0 

0.0 

0.0 

0.0 

B.2  SHEL  EOS  MODELS  GENERIC2  (GRAN1 C  AT  4,500  M/S)  AND 
GENERICS  (GRAN1E  AT  4,500  M/S). 


123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789- 


GRANIC  1  1 

GRANIA  with 

new  FAN 

2.624 

2.650 

0.04 

0.0 

~120.E-06 

0.30 

0.7 

.03 

0.0 

1.000 

0.0 

1.00 

0.0 

0.020 

0.09 

0.0 

0.0 

0.0 

0.0 

0.05 

1.0 

3.0 

0.00 

0.0 

4.30 

.007 

.13 

2.0 

2.3 

0.25 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0125 

0.0001 

0.3 

0.0 

0.0 

0.0 

0.0 

0.0 
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IE  1  1 

GRANIA  with 

new  FAN 

and  EXP 

yield  surface 

2.624 

2.650 

0.04 

0.0  • 

-120.E-06 

0.30 

0.7 

.03 

0.0 

1.000 

0.0 

1.00 

0.0 

0.020 

0.09 

0.0 

0.0 

0.0 

0.0 

0.05 

1.0 

3.0 

0.00 

0.0 

4.30 

.007 

.13 

2.0 

2.3 

0.25 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0125 

0.0001 

30.0 

0.0 

0.0 

0.0 

0.0 

0.0 

B.3  JWL  HIGH  EXPLOSIVE  MODEL  FOR  QM-100/R  (LLNL  FIT). 

123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789- 
QM10R6  4  1  2ND  CUT  LLNL  JWL  FOR  QMIOOR  8/2/93 

1.509  1.509  0.043075  0.043075  -1.0 

0.56000  11.025251  6.50  0.744978  2.8350  0.74200  .0650  0.20 

The  JWL  constants  of  the  above  line  are  as  follows: 

w  A(Mb)  R1  B(Mb)  R2  DetV(cni/ms)  Eo(Te/cc)  Pcj  (Mb) 


B.4  MCIST  EOS  MODEL  CON5K4  FOR  HJ7C  CONCRETE. 


123456789-123456789-123456789-123456789-123456789-123456789-123456789-123456789- 


CON5K4  5 

1  NEW  5KSI 

CONCRETE 

FIT  TO  MN 

CAL  DATA 

3/19/93  VML= 

5KB 

2.279 

2.279 

0.03 

0.0  - 

-24.0E-06 

7.348 

0.363 

0.0143028 

0.206272 

.400 

.180 

O.lE-7 

1.0 

0.3  .35 

0.0 

0.0 

1.0 

0.0 

0.0 

0.0 

0.0  24. 

3.0 

0.03 

0.0 

0.0 

0.0 

0.0 

0.0  0.95 

0.26 

0.18 

1.0 

0.0 

0.0 

0.0 

0.0  2.89E-3 

2.89E-3 

49.5E-6 

1.484 

0.0 

0.0 

0.0 

2.0E--6  1.484 
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CONC5K4  MClST  DATA  FILE  (REQUIRED  FOR  MOIST  MODEL) 


O.OOOOOE+OO, 

1.08696E-04, 

2.17392E-04, 

1.08702E-03, 

1.28700E-03, 

1.48700E-03, 

1.68700E-03, 

1.88700E-'03, 

2.08700E-03, 

2.28700E-03, 

2.48700E-03, 

2.68700E-03, 

2.88700E-03, 

3.08700E-03, 

3.28700E-03, 

3.4.8700E-03, 

3.68700E-03, 

3.88700E-03, 

4.08700E-03, 

4.28700E-03, 

4.48700E-03, 

4.68700E-03, 

4.88700E-03, 

5.08700E-03, 

5.28700E-03, 

5.48700E~03, 

5.68700E-03, 

5.88700E-03, 

6.08700E~03, 

6.28700E-03, 

6.48700E-03, 

6.68700E-03, 

6.88700E-03, 

7.08700E~03, 

7.28700E-03, 

7.48700E-03, 

7.68700E-03, 

7.88700E-03, 

8.08700E-03, 


O.OOOOOE+00 

2.83784E-02 

5.67568E-02 

2.83800E-01 

3,23700E-01 

3.57400E-01 

3.86300E-01 

4.12200E-01 

4.36300E-01 

4.60300E-01 

4.84500E-01 

5.08600E-01 

5.32300E-01 

5.55100E-01 

5.76700E-01 

5.96800E-01 

6.15300E-01 

6.32500E-01 

6.48800E-01 

6.64500E-01 

6.79900E-01 

6.95300E-01 

7.11100E-01 

7,27000E-01 

7.43000E-01 

7.58800E~01 

7.74300E-01 

7.89300E-01 

8.03700E-01 

8.17500E-01 

8.30700E-01 

8.43500E-01 

8.55800E-01 

8.67800E-01 

8.79500E-01 

8.91000E-01 

9.02300E-01 

9.13600E-01 

9.24700E-01 


8.28700E-03 

8.48700E-03 

8.68700E-03 

8.88700E-03 

9.08700E-03 

9.28700E-03 

9.48700E-03 

9.68700E-03 

9.88700E-03 

1.00900E-02 

1.02900E-02 

1.04900E-02 

1.06900E-02 

1.08900E-'02 

1.10900E-02 

1.12900E-02 

1.14900E-02 

1.16900E-02 

1.18900E-02 

1.20900E-02 

l,22900E-02 

1.24900E-02 

1.26900E-02 

1.28900E-02 

1.30900E-02 

1.32900E~02 

1.34900E-02 

1.36900E-02 

1.38900E-02 

1.40900E-02 

1.42900E-02 

1.44900E-02 

1.46900E-02 

1.48900E-02 

1.50900E-02 

1.52900E-02 

1.54900E-02 

1.56900E-02 

1.58900E-02 


9.35700E-01 

9.46600E-01 

9.57400E-01 

9.68000E-01 

9.78500E-01 

9.88800E-01 

9.98900E-01 

1.00900E+00 

1.01900E+00 

1.02800E+00 

1.03800E+00 

1.04700E+00 

1.05700E+00 

1.06600E+00 

1.07500E+00 

1.08400E+00 

1.09300E+00 

1.10200E+00 

l,11100E+00 

1.12000E+00 

1.12900E+00 

1.13800E+00 

1.14700E+00 

1.15600E+00 

1.16400E+00 

1,17300E+00 

1.18106e+00 

1.19000E+00 

1.19800E+00 

1.20600E+00 

1.21400E+00 

1.22200E+00 

1.22900E+00 

1.23700E+00 

1.24500E+00 

1.25200E+00 

1.26000E+00 

1.26700E+00 

1.27500E+00 
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1.60900E-02, 

1.28300E+00 

2.42900E-02, 

1.61000E+00 

1.62900E-02, 

1.29000E+00 

2.44900E-02, 

1.61800E+00 

1.64900E-02, 

1.29800E+00 

2.46900E-02, 

1.62500E+00 

1.66900E-02, 

1.30600E+00 

2.48900E-02, 

1.63300E+00 

1.68900E-02, 

1.31400E+00 

2.50900E-02, 

1.64000E+00 

1.70900E-02, 

1.32200E+00 

2.52900E-02, 

1.64800E+00 

1.72900E-02, 

1.32900E+00 

2.54900E-02, 

1.65500E+00 

1.74900E-02, 

1.33700E+00 

2.56900E-02, 

l,66300E+00 

1.76900E-02, 

1.34500E+00 

2.58900E-02, 

1.67000E+00 

1.78900E-02, 

1.35300E+00 

2.60900E-02, 

1.67800E+00 

1.80900E-02, 

1.36100E+00 

2.62900E-02, 

1.68500E+00 

1.82900E-02, 

1.36900E+00 

2.64900E-02, 

l,69200E+00 

1.84900E-02, 

1.37800E+00 

2.66900E-02, 

1.70000E+00 

1.86900E-02, 

1.38600E+00 

2.68900E-02, 

1.70700E+00 

1.88900E-02, 

1.39400E+00 

2.70900E-02, 

1.71500E+00 

1.90900E-02, 

1.40200E+00 

2.72900E-02, 

1.72200E+00 

1.92900E-02, 

1.41000E-I-00 

2.74900E-02, 

1.73000E+00 

1.94900E-02, 

1.41800E+00 

2.76900E-02, 

1.73700E+00 

1.96900E-02, 

1.42600E+00 

2.78900E-02, 

1.74500E+00 

1.98900E-02, 

1.43500E+00 

2.80900E-02, 

1.75200E+00 

2.00900E-02, 

1.44300E+00 

2.82900E-02, 

1.76000E+00 

2.02900E-02, 

1.45100E+00 

2.84900E-02, 

1.76700E+00 

2.04900E-02, 

1.45900E+00 

2.86900E-02, 

1.77500E+00 

2.06900E-02, 

1.46700E+00 

2.88900E-02, 

1.78300E+00 

2.08900E-02, 

1.47600E+00 

2.90900E-02, 

1.79000E+00 

2.10900E-02, 

1.48400E+00 

2.92900E-02, 

1.79800E+00 

2.12900E-02, 

1.49200E+00 

2.94900E-02, 

1.80600E+00 

2.14900E-02, 

1.50000E+00 

2.96900E-02, 

1.81400E+00 

2.16900E-02, 

1.50800E+00 

2.98900E-02, 

1.82200E+00 

2.18900E-02, 

1.51600E+00 

3.00900E-02, 

1.83000E+00 

2.20900E-02, 

l,52400E+00 

3.02900E-02, 

1.83700E+00 

2.22900E-02, 

-1.53200E+00 

3..04900E-02, 

1.84500E+00 

2.24900E-02, 

1.54000E+00 

3.06900E-02, 

1.85300E+00 

2,26900E-02, 

l,54800E+00 

3 .08900E-02, 

1.86100E+00 

2.28900E-02, 

1.55600E+00 

3.10900E-02, 

1.86900E+00 

2.30900E-02, 

1.56400E+00 

3 .12900E-02, 

1.87700E+00 

2.32900E-02, 

1.57200E+00 

3.14900E-02, 

1.88500E+00 

2.34900E-02, 

1.58000E+00 

3 .16900E-02, 

1.89300E+00 

2.36900E-02, 

1.58700E+00 

3.18900E-02, 

1.90100E+00 

2.38900E-02, 

1.59500E+00 

3.20900E-02, 

1.90900E+00 

2.40900E-02, 

1.60300E+00 

3.22900E-02, 

1.91700E+00 
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3,24900E-02, 

3.26900E-02, 

3.28900E-02, 

3.30900E-02, 

3.32900E-02, 

3.34900E-02, 

3.36900E-02, 

3.38900E-02, 

3.40900E-02, 

3.42900E-02, 

3.44900E-02, 

3.46900E-02, 

3.48900E-02, 

3.50900E-02, 

3.52900E-02, 

3.54900E-02, 

3.56900E*-02, 

3.58900E-02, 

3.60900E-02, 

3.62900E-02, 

3.64900E-02, 

3.66900E-02, 

3.68900E-02, 

3.70900E-02, 

3.72900E-02, 

3.74900E>02, 

3 .76900E-02, 

3..78900E-02, 

3.80900E-02, 

3.82900E-02, 

3.84900E-02, 

3 .86900E-02, 

3 .88900E-02, 

3.90900E-02, 

3.92900E~02, 

3 .94900E-02, 

3.96900E-02, 

3.98900E-02, 

4.00900E-02, 

4.02900E-02, 

4.04900E-02, 


1.92400E+00 

1.93200E+00 

1.94000E+00 

1.94700E+00 

1.95500E+00 

1.96300E+00 

1.97000E+00 

1.97700E+00 

1.98500E+00 

1.99200E+00 

1.99900E+00 

2.00700E+00 

2.01400E+00 

2.02100E+00 

2.02800E+00 

2.0350QE+00 

2.04200E+00 

2.04900E+00 

2.05600E+00 

2.06300E+00 

2.07100E+00 

2.07800E+00 

2.08500E+00 

2 .09200E+00 

2.09900E+00 

2.10700E+00 

2.11400E+00 

2.12200E+00 

2.12900E+00 

2.13700E+00 

2.14400E+00 

2.15200E+00 

2.16000E+00 

2.16800E+00 

2.17600E+00 

2.18400E+00 

2.19200E+00 

2.20000E+00 

2.20900E+00 

2.21700E+00 

2.22600E+00 


4.06900E-02 

4.08900E-02 

4.10900E-02 

4.12900E-02 

4.14900E-02 

4.16900E-02 

4.18900E-02 

4.20900E-02 

4.22900E-02 

4.24900E-02 

4,26900E-02 

4.28900E-02 

4.30900E-02 

4.32900E-02 

4.34900E-02 

4.36900E-02 

4.38900E-02 

4.40900E-02 

4.42900E-02 

4.44900E-02 

4.46900E-02 

4.48900E-02 

4.50900E-02 

4,52900E-02 

4.54900E-02 

4.56900E-02 

4.58900E-02 

4.60900E-02 

4.62900E-02 

4.64900E-02 

4.66900E-02 

4.68900E-02 

4.70900E-02 

4.72900E-02 

4.74900E-02 

4.76900E-02 

4.78900E-02 

4.80900E-02 

4.82900E-02 

4.84900E*>02 

4.86900E~02 


2.23400E+00 

2.24300E+00 

2.25200E+00 

2.26100E+00 

2.27000E+00 

2.27900E+00 

2.28800E+00 

2.29700E+00 

2.30600E+00 

2.31500E+00 

2.32400E+00 

2.33300E+00 

2.34300E+00 

2.35200E+00 

2.36100E+00 

2,37000E+00 

2.37900E+00 

2.38800E+00 

2.39700E+00 

2.40600E+00 

2.41500E+00 

2,42400E+00 

2.43300E+00 

2,44200E+00 

2.45000E+00 

2.45900E+00 

2.46800E+00 

2.47700E+00 

2.48600E+00 

2.49500E+00 

2.50400E+00 

2.51200E+00 

2.52100E+00 

2.53000E+00 

2.53900E+00 

2.54800E+00 

2.55700E+00 

2.56600E+00 

2.57500E+00 

2.58400E+00 

2.59300E+00 
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4.88900E-02, 

2.60200E+00 

5.70900E-02, 

2.99300E+00 

4.90900E-02, 

2.61100E+00 

5.72900E-02, 

3.00200E+00 

4.92900E-02, 

2 .62000E+00 

5.74900E-02, 

3.01100E+00 

4.94900E-02, 

2.63000E+00 

5.76900E-02, 

3.02100E+00 

4.96900E-02, 

2.63900E+00 

5.78900E-02, 

3.03000E+00 

4.98900E-02, 

2.64800E+00 

5.80900E-02, 

3.03900E+00 

5.00900E’-02, 

2.65700E+00 

5.82900E-02, 

3.04900E+00 

5.02900E-02, 

2.66700E+00 

5.84900E-02, 

3.05800E+00 

5.04900E-02, 

2.67600E+00 

5.86900E-02, 

3.06700E+00 

5.06900E-02, 

2.68500E+00 

5.88900E-02, 

3.07600E+00 

5.08900E-02, 

2.69500E-H00 

5.90900E-02, 

3.08600E+00 

5.10900E-02, 

2.70400E+00 

5.92900E-02, 

3.09500E+00 

5.12900E-02, 

2.71400E+00 

5.94900E-02, 

3.10400E+00 

5.14900E-02, 

2.72300E+00 

5.96900E-02, 

3.11400E+00 

5.16900E-02, 

2.73300E+00 

5.98900E-02, 

3.12300E+.00 

5.18900E-02, 

2.74200E+00 

6.00900E-02, 

3.13200E+00 

5.20900E-02, 

2.75200E+00 

6.02900E-02, 

3.14200E+00 

5.22900E-02, 

2.76100E+00 

6.04900E-02, 

3.15100E+00 

5.24900E-02, 

2.77100E+00 

6.06900E-02, 

3 .16000E+00 

5.26900E-'02, 

2.78100E+00 

6.08900E-02, 

3.17000E+00 

5.28900E-02, 

2.79000E+00 

6.10900E-02, 

3.17900E+00 

5.30900E-02, 

2,80000E+00 

6.12900E-02, 

3.18900E+00 

5.32900E-02, 

2.81000E+00 

6.14900E-02, 

3.19800E+00 

5.34900E-02, 

2.82000E+00 

6.16900E-02, 

3.20800E+00 

5.36900E-02, 

2.82900E+00 

6.18900E-02, 

3.21700E+00 

5.38900E-02, 

2,83900E+00 

6.20900E-02, 

3.22700E+00 

5.40900E-02, 

2.84900E+00 

6.22900E-02, 

3.23600E+00 

5.42900E-02, 

2.85900E+00 

6.24900E-02, 

3.24600E+00 

5.44900E-02, 

2.86800E+00 

6.26900E-02, 

3.25600E+00 

5.46900E-02, 

2.87800E+00 

6.28900E-02, 

3.26500E+00 

5.48900E-02, 

2.88800E+00 

6.30900E-02, 

3.27500E+00 

5.50900E-02, 

2 .89700E+00 

6.32900E-02, 

3.28500E+00 

5.52900E-02, 

2.90700E+00 

6.34900E-02, 

3.29500E+00 

5.54900E-02, 

2.91700E+00 

6.36900E-02, 

3.30500E+00 

5.56900E-02, 

2.92600E+00 

6.38900E-02, 
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